Phase diagrams of a classical two-dimensional Heisenberg antiferromagnet with 

single-ion anisotropy 
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A classical variant of the two-dimensional anisotropic Heisenberg model reproducing inelastic 
neutron scattering experiments on La5CagCu2404i [M. Matsuda et at, Phys. Rev. B 68, 060406(R) 
(2003)] is analysed using mostly Monte Carlo techniques. Phase diagrams with external fields 
parallel and perpendicular to the easy axis of the anisotropic interactions are determined, including 
antiferromagnetic and spin-flop phases. Mobile spinless defects, or holes, are found to form stripes 
which bunch, debunch and break up at a phase transition. A parallel field can lead to a spin-flop 
phase. 
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I. INTRODUCTION 

The compounds (La,Ca)i4Cu2404i have been 
studied experimentally rather extensively in recent 
years >i»2i2diSiSiL2iS They display interesting low-di- 
mensional magnetic properties arising from CU2O3 
two- leg ladders and CuC>2 chains. In addition, for 
Lai4_ x Ca Z :Cu2404i there is an intriguing interplay 
between spin and charge ordering due to hole doping 
when x > 8. 

In particular, La 5 CagCu2404i exhibits, at low tem- 
peratures and small fields, antiferromagnetic long range 
order associated with the CuC>2 chains which are ori- 
ented along the c axis. The Cu 2+ ions in those chains 
carry a spin-1/2. The spins are ordered ferromagnetically 
in the chains, while the interchain coupling in the ac— 
planes is antiferromagnetic. The magnetic interactions 
between the ac-planes are believed to be very weak. The 
couplings in the ac-planes are anisotropic with an easy 
axis along the b axis, i.e. the magnet has an Ising-likc 
characteri5i£i2i& Holes may originate from the La and Ca 
ions, transforming Cu ions into spinless quantities, with 
a hole content of about 10 percent^ 

Experiments on La5CagCu2404i include thermody- 
namic measurements on the specific heat, magnetization, 
and susceptibility^ as well as electron spin resonance 6 
and neutron scattering^^ Motivated by the experimen- 
tal findings, different models have been proposed and 
studied. Firstly, a two-dimensional Ising model has been 
introducedfi^i^ where the spins correspond to the Cu 2+ 
ions, and mobile spinless defects mimic the holes. In ad- 
dition to the ferromagnetic intrachain and antiferromag- 
netic interchain couplings between neighboring spins, a 
rather strong and antiferromagnetic exchange between 
next-nearest chain spins separated by a hole is assumed. 
The model has been shown to describe, at low tem- 
peratures, antiferromagnetic domains separated by quite 
straight defect lines which break up at a phase transi- 
tion where also the long range magnetic order gets de- 
stroyed. The stripe instability is caused by an effectively 
attractive interaction between the defects mediated by 



the antiferromagnetic interchain couplings. 

Even more recently, a two-dimensional anisotropic 
Heisenberg model has been shown by Matsuda et al. 
to reproduce the measured spin-wave dispersions which 
supposedly result from the collective spin excitations of 
the Cu 2+ ions in the ac-planes£ Our subsequent Monte 
Carlo simulations on the classical variant of the model 
seem to indicate, however, that some of its thermo- 
dynamic properties tend to deviate from experimental 
findings^ In particular, in an external field along the b 
axis, at low temperatures, the field dependence of the sus- 
ceptibility of the anisotropic Heisenberg model disagrees 
with the measured behavior. The disagreement can prob- 
ably not be resolved by invoking quantum effects. In par- 
ticular, the critical temperature of the model seems to be 
significantly lower than the measured one, as can be seen 
by comparing results on related quantum and classical 
models(ifl'ii*i£ Indeed, in this paper we will show that 
the qualitative properties of the phase diagram of the 
model are not affected qualitatively by quantum effects. 

In any event, the classical variant of the model of 
Matsuda et al. deserves to be studied in more detail. 
Apart from the previous qualitative comparison^ with 
experimental data, the classical model is of genuine the- 
oretical interest as well. Perhaps most interestingly, the 
phase diagram in the (temperature,field)-plane is worth 
studying in the context of two-dimensional anisotropic 
Heisenberg models. Due to the anisotropy, there is a 
spin-flop phase when applying a sufficiently high ex- 
ternal field along the b axis. In two dimensions, that 
phase has interesting properties as it is believed to be 
of Kosterlitz-Thouless type with spin correlations de- 
caying algebraically with distance. Furthermore, the 
boundary line of the antiferromagnetic phase as well as 
the transition between the antiferromagnetic and the 
spin-flop phase has been discussed controversially for 
anisotropic nearest-neighbor Heisenberg models. In three 
and higher dimensions, there is a bicritical point where 
the antiferromagnetic, the spin-flop, and the paramag- 
netic phases meet4& As it is well known, such a point 
could occur in two dimensions only at zero tempera- 
ture (Mermin- Wagner theorem). Actually, different see- 
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narios have been proposed^ 7 - However, simulations and 
other analyses on classical models lead to controversial 
results iA 7 *^ Fairly recently, a S = 1/2 quantum version 
has been simulated suggesting a topology of the phase 
diagram with a tricritical and a critical end point. 19 Our 
work will deal with this aspect. Note that quasi two- 
dimensional anisotropic Heisenbcrg antifcrromagnets ex- 
hibiting a spin-flop phase have attracted much exper- 
imental attention already some years agOi22i2ii2£ Well- 
known examples are Ptf^MnCLt and K2M11F4. The ap- 
proach of our study may be also useful for the correct 
theoretical analysis of models describing these materials. 

Our main emphasis will be on a classical variant of the 
model by Matsuda et al. obtained from the spin- wave dis- 
persions of La5CagCu2404i. In addition, in an attempt 
to include possible effects due to the holes, we shall ex- 
tend the model by allowing for mobile defects following 
the previous considerations^ In fact, the defects again 
form stripes which are destabilized at a phase transition. 
Bunching and debunching of the stripes are novel fea- 
tures. Effects of an external field along the b axis will 
also be considered. 

The layout of the article is as follows: In the next sec- 
tion, the model will be introduced. Its phase diagrams, 
without defects and applying external fields parallel and 
perpendicular to the easy axis of the magnetic interac- 
tions, will be presented in Sec. lIIII The possible influence 
of defects on thermal properties will be discussed in sec- 
tion ^3 followed by a summary. 

II. THE MODEL, SIMULATIONS, AND 
QUANTITIES OF INTEREST 

Following Matsuda et al., 9 the magnetic properties 
of La5CagCu2404i depend on the Cu 2+ ions located in 



where S; >m = (Sf m , Sf m , Sf m ) denotes the spin at the l- 
th site in the m-th chain; a refers to differently oriented 
magnetic fields, with a = z for the field parallel to the 
easy axis, and a — x or y for a field perpendicular to it, 
i.e. along the a or c axis. 

Extending this Hamiltonian of Matsuda et al, we in- 
troduce defects S;, m = 0, induced by the holes originating 
from the La and Ca ions m& We neglect direct interactions 
between the defects, and there are no couplings between 
a defect and a spin. Next-nearest neighbor spins in the 



the ac-planes, having a centered rectangular geometry 
as depicted in Fig. ^ Based on their spin-wave anal- 
ysis, the spins (S = 1/2) of the ions couple along the 
Cu02 chains, i.e. along the c axis, through nearest neigh- 
bor, J c i, and next-nearest neighbor, J C 2, exchange con- 
stants, with J c i = —0.2 meV being antiferromagnetic and 
J C 2 = 0.18 meV being ferromagnetic. The ferromagnetic 
ordering in the chains is due to the strong antiferromag- 
netic interchain couplings, see also Fig.^ Jaci = —0.681 
meV refers to the two nearest neighbors in the adjacent 
chain, and J aC 2 = 0.5J ac i = —0.3405 meV denotes the 
couplings to the two next-nearest neighbors. 



Importantly, there is an uniaxial exchange anisotropy 
favoring alignment of the spins along the b axis. Its con- 
tribution to the different couplings cannot be determined 
in the spin-wave analysis, and only its integral effect on 
the gap in the dispersion is quantified. 9 When going over 
to a classical description with spins of fixed length, say, 
one, the total anisotropy may be mimicked by a single- 
ion term. Such a term would be, of course, merely a triv- 
ial constant and unphysical for a quantum system with 
S= 1/2. The single-ion term, avoiding ambiguities in dis- 
tributing the anisotropy among the different couplings, is 
quite reasonable in the classical variant. It leads to about 
the same transition temperature without field as various 
exchange anisotropies of the same overall magnitude. 15 
In addition, the entire phase diagram seems to be not 
affected by the type of anisotropy. The anisotropy has 
been estimated from the gap energy to be D = —0.211 
meV, 9 Thence the Hamiltonian can be written a o 9 ' 1 ' 5 " 



(1) 



same chain, Si m and Si±2 iTn , with a defect in between 
them are coupled by the exchange constant Jo (replac- 
ing J C 2), which we presume to be antiferromagnetic and 
rather strong as before, i.e. of several meV>i^ Specifically, 
we choose Jo = —6.25 meV. Thence, the Hamiltonian JJJ 
is augmented by the term 

Ti-d = — Jo ^ S^„ l S'; + 2, m (l — Hi+l, m ), (2) 

(i,m) 

with ni^m = 0, 1 being the occupation variable of a spin 



ti — —Jci Sl,mSl+l,m — Jc2 S/,mS;+2,m — Jael &l,m (Sl,m+1 + S; + i, m +i) 

(l,m) (,h m ) 

- Jac2 ^ S; )m (S;_i, m+ i + S;+2, m +l) + D } J {Si m ) 2 - H a ^ S" m , 



(l,m) 



(l,m) 



(l,m) 



r 



3 




(Z + l,m + l) 
(l + 2,m) 



FIG. 1: Sketch of the the magnetic interactions for the two- 
dimensional anisotropic Heisenberg model of Matsuda et aim 
The dots denote the sites of the Cu 2+ ions in the ac-plane of 

Las Cag CU24 O41 • 



at site (l,m). The defects are allowed to hop to a neigh- 
boring site in a chain, transforming the spin at that site 
into a defect and leaving a new spin with arbitrary orien- 
tation at its initial site. The probability of such a process 
is determined by the Boltzmann factor of the change in 
the magnetic energy, Eqs. and (J2J), associated with 
the hop»i£ It is easily incorporated in Monte Carlo sim- 
ulations. 

As before^ we assume that defects are separated 
along the chain by at least one spin. The number of de- 
fects in each chain will be taken to be 10 percent of the 
number of sites in that chain. The defect concentration 
is then close to that in La5CagCu2404i AiS 

We shall study the model, with and without defects, 
using, apart from ground state calculations, standard 
Monte Carlo techniques with local elementary processes 
changing at randomly chosen sites the spin orientation 
and moving defects to neighboring sites. In the sim- 



ulations, we consider lattices with the same number L 
of chains and of sites per chain, i.e. with a total of L 2 
sites. Full periodic boundary conditions are employed. 
To study finite-size effects allowing extrapolations to the 
thermodynamic limit, L — > 00, we consider typically 
sizes ranging from L=10 to L=200. Each run consists of 
at least 10 6 and up to 5 x 10 6 Monte Carlo steps per site. 
To obtain averages and error bars, we take into account 
the results of up to 10 realizations using different random 
numbers. 

We compute both quantities of direct experimental in- 
terest as well as other quantities which enable us to de- 
termine critical properties and the phase transition lines. 
In particular, we recorded the specific heat C, both from 
the fluctuations and from the temperature derivative of 
the energy per site E. In the absence of defects, various 
magnetizations were computed. Especially, we recorded 

(i) the a-component of the magnetization to study the 
response to a field in the a direction, with a = x,y, z, 



< M a >- 



~< E ^"n 

(l,m) 



>/L 2 



(3) 



(ii) the z-component of the absolute value of the stag- 
gered magnetization M z s and the square of the staggered 
magnetization to describe the order in the antiferromag- 
netic phase, 



< IMf I >=< 



E 

(all /, m even) 



(all I, m odd) 



(4) 

(summing separately over sites in even and odd chains), 
and similarly for < (M|) 2 >; and 

(iii) the square of the transverse sublattice magnetization 
to describe the Kosterlitz-Thouless character of the spin- 
flop phase when applying a field along the easy axis, 



J 



< M. 



■HJ 



> = < 



E 



[( E 



SI 



l.rn 



(all I, m even) 



E * 



(all /, m odd) 



I, in 



>/(L 4 /2). 



(5) 



In addition, we recorded the magnetic susceptibilities 
X a , which may be computed from the fluctuations or 
field derivatives of the corresponding magnetizations, 
< M a >, and the (finite lattice) staggered susceptibil- 
ity Xs defined by 



X* = L 2 (<(M*) 2 >-<|AC|> 2 )/2. 



(0) 



To identify the type of transition from the antiferromag- 
netic to the paramagnetic and the spin-flop phases, the 
fourth-order, size dependent cumulant of the staggered 
magnetization, the Binder cumulant 23 Ul, is rather use- 



ful: 



Ur 



1- < {M z s f > /(3 < {M z s f > 2 ) 



(7) 



where < (M|) 4 > is defined in analogy to < (M z ) 2 >. 

In the presence of defects we studied, apart from the 
specific heat, magnetization and susceptibility, also mi- 
croscopic quantities describing the topology and stability 
of the defect stripesfi 3 ^ including the average minimal 
distance between defects in adjacent chains and the den- 
sity of defect pairs. Results will be discussed in Sec. IIVI 
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III. PHASE DIAGRAMS IN THE ABSENCE OF 
DEFECTS 

We analysed the anisotropic Heisenberg model of Mat- 
suda et al& applying external fields along the easy axis 
(H z ) and perpendicular to it (say, H x ), and varying the 
temperature, see the Hamiltonian 0). The resulting 
phase diagrams are depicted in Fig. |21 



g- -o field along easy axis (IT) 
h— a field perp. to easy axis (H ) 



paramagnetic 




paramagnetic 



antiferromagnetic 



0.74 0.76 



0.78 



k B T 



0.8 



0.82 0.84 



FIG. 2: Phase diagram of the model without defects, with 
fields parallel, H z , and perpendicular, H x , to the easy axis, 
(a) Global phase diagram, (b) details in the (T,H Z ) plane. 
Crosses denote approximate locations of the tricritical point 
and the critical end point (see text). 

In the case of a field H x > perpendicular to the easy 
axis one encounters, at zero temperature and small fields, 
H x < H% m , an antiferromagnetic ground state with the 
non-zero, field-dependent z-component of the spins in 
each chain pointing in the same direction and alternat- 
ing sign from chain to chain, Sf m = —S' m+1 . At zero 
temperature and larger fields, H x > H£ m , the magnetic 
field term dominates, and the spins are aligned along the 
direction of the field, Sf m =l. The critical field H?™ is 
readily calculated. Inserting the values of the intrachain, 
J c i and J C 2, and interchain coupling constants, J ac \ and 



Jaci-, as well as the spin anisotropy, D, as stated in the 
preceeding section, one gets H? m — 8.594.. meV. Cer- 
tainly, this is an artificial unit, which had to be tran- 
scribed into the standard unit Tesla, taking into account 
the ^-factor and the actual spin value, when comparing 
results for non-vanishing fields quantitatively to experi- 
mental findings. However, in the context of our analysis 
the artificial unit will be sufficient. In the following, the 
unit "meV" will be suppressed in all expressions for the 
energy, ks times temperature (fcgT), and the magnetic 
field. 

At non-zero temperatures, a critical line arises from 
(T = 0,H X = HP™) ending at (T = T N ,H X = 0), see 
Fig. The transition separates the ordered antiferro- 
magnetic phase with a non-zero staggered magnetization 
< \M"\ >, see Eq. (0J, from the disordered (paramag- 
netic) phase, where < \M"\ >= 0. The phase tran- 
sition is expected to be continuous and of Ising type, 
i.e. with the well-known critical exponents of the two- 
dimensional Ising model. The critical line has been ob- 
tained by fixing either the temperature and varying the 
field or by fixing the field and varying the temperature. 
Then standard finite-size analyses on the peak positions 
of the specific heat, T%(L), were done. Indeed, these 
positions approach, for sufficiently large system sizes, 
the critical temperature T c of the infinite system like 
T C — T^(L) oc 1/L, which is consistent with the transition 
belonging to the Ising universality class. For illustrative 
purposes, some raw data on the specific heat are shown 
in Fig. |21 The thermal behavior of the staggered magne- 
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FIG. 3: Specific heat vs. temperature at fixed fields H x — 
H z = (solid lines), H x = 4.0 (dotted lines), and H z = 4.0 
(dashed lines), for systems, without defects, of various sizes 
L. 

tization is shown, for a few selected examples, in Fig. 0] 

In the case of an external field H z > along the easy 
axis one obtains a more complex and more interesting 
phase diagram, see Fig. [3 In the ground state (T = 0), 



5 




0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 
k B T 

FIG. 4: Staggered magnetization vs. temperature for various 
fields parallel (H z ) and perpendicular {H x ) to the easy axis, 
simulating systems, without defects, of size L — 100. 

one has to distinguish two critical fields, Hf and Hf m . 
For H z < Hf, the antiferromagnetic structure, as de- 
scribed above, has the lowest energy. At larger fields, 
Hf < H z < -H]? 111 , the spin-flop state is stable. There 
the z-component of the spins in all chains acquires the 
same, field-dependent value S Z (H Z ) > 0. The planar, 
sy-components of the spins are aligned parallel to each 
other in each chain, pointing in an arbitrary direction 
due to the rotational invariance of the interactions in the 
xy-plane, see Eq. Q). The iry-components of spins in 
neighboring chains point in the opposite direction be- 
cause of the antiferromagnetic interchain couplings. At 
H > Hf m , one has a ferromagnetic ordering with S z = 1. 

For the set of couplings obtained from the spin-wave 
analysis, the critical fields are Hf = 1.808.. and Hf m = 
7.75. At H z = Hf, the ^-component takes the value 
S z (Hf) = 0.233.., corresponding to an angle of 76.5.. 
degrees formed by the z axis and the orientation of the 
spins. 

The complete phase diagram in the (T, H z ) plane con- 
sists of the antiferromagnetic, the spin-flop and the para- 
magnetic or disordered states, see Figs. [2^ and|2p. The 
antiferromagnetic phase exhibits long-range order with 
the staggered magnetization M z as order parameter. 
The spin-flop phase has been argued to be of Kosterlitz- 
Thouless character jiSiiLiSiiS where transverse spin corre- 
lations, i.e. < Sf m Sf, m , +Sf m Sf, m , >, decay algebraically 
with distance — I') 2 + (m — m') 2 . Accordingly, the 
transverse sublattice magnetization, see Eq. (J5J, being 
the order parameter of the spin-flop phase in three di- 
mensions, is expected to behave for T > and sufficiently 
large systems as 

< Ml >cx L-« (8) 

with g approaching 1/4 at the transition from the spin- 
flop to the paramagnetic phased and g — 2 in the para- 
magnetic phase. Thence, the order parameter vanishes 



in the Kosterlitz-Thouless phase as L — ► oo at all tem- 
peratures T > 0. Of course, in the disordered phase spin 
correlations decay exponentially with distance. 

While the existence of these phases for weakly an- 
isotropic Heisenberg antiferromagnets in two dimen- 
sions is undisputed, basic aspects of the topology of 
the phase diagram and especially the transitions be- 
tween the antiferromagnetic phase and the spin-flop as 
well as the paramagnetic phases have been discussed 
controversiallyjiZii2ii£ and they may, indeed, depend on 
details of the model. 

We determined the boundary line of the antiferromag- 
netic phase by monitoring especially the specific heat C, 
the (square of the) staggered magnetization, < \M Z \ > 
and < (M z ) 2 >, the staggered susceptibility, \ z s , an d 
the Binder cumulant, Ul, Eq. @. A few raw data for 
the specific heat and the staggered magnetization are in- 
cluded in Figs. and 

The transition from the antiferromagnetic to the disor- 
dered phase at low fields and high temperatures is contin- 
uous and of Ising type. Its location, as displayed in Figs. 
|2t and|2p, follows from finite-size analyses of the various 
physical quantities. The data are consistent with a loga- 
rithmic divergence of the specific heat as well as with the 
well-known Ising values for the critical exponents of the 
order parameter, (3 — 1/8, and of the staggered suscepti- 
bility, 7 = 7/4. 

More interestingly, the transition from the antifer- 
romagnetic to the paramagnetic phase eventually be- 
comes first order when increasing the field and lower- 
ing the transition temperature, with a tricritical point 
at k B T tI = 0.79 ± 0.015 and Hf = 1.665 ± 0.01. The 
boundary line of the antiferromagnetic phase remains 
to be first order at lower temperatures when separat- 
ing the antiferromagnetic and the spin-flop phase. The 
Kosterlitz-Thouless line separating the spin-flop phase 
from the paramagnetic (disordered) state hits the bound- 
ary of the antiferromagnetic phase in a critical end point 
at k B T ccp = 0.75 ± 0.015 and H™ p = 1.675 ± 0.01, 
see Fig. [2J3. Note that the phase diagram has qual- 
itatively the same topology as the one suggested for 
the spin-1/2 quantum version of the standard nearest- 
neighbor antiferromagnet with exchange anisotropy in 
two dimensions^ in agreement with the classical ver- 
sion of that model as welliS Therefore, we conclude that 
quantum effects are of minor importance for the main 
features of the phase diagram. 

The tricritical point may be located by studying the 
Binder cumulant. In the thermodynamic limit the value 
of the cumulant at the transition point, C/l =00 , is known 
to depend on the type and universality class of the tran- 
sition. In simulations, Ul=<x> can be estimated from sys- 
tematic finite-size extrapolations of the intersection val- 
ues of the Binder cumulant for different system sizes L\ 
and L 2 , U Ll = Ul 2 = U(Li,L 2 ). 23 In Fig. we depict 
results for [7(80, 100), obtained usually at fixed tempera- 
ture and varying the field in the vicinity of the boundary 
of the antiferromagnetic phase. Obviously, J7(80, 100) is 
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nearly constant at high temperatures, U ~ 0.6, with a 
fairly rapid change around k^T w 0.80. This finding and 
further finite-size analyses on U(L\,L?) for other sys- 
tem sizes allow us, indeed, to approximately locate the 
tricritical point which separates the transition of Ising 
type, where Ul=oo = 0.6106. .^21 and the transition of 
first order. Note that the value of Ul=oo may be slightly 
affected due to the interactions J c i, J C 2 and J aC 2- If 
only the predominant coupling J ac \ were non-zero, the 
model is easily seen to be equivalent to a nearest-neighbor 
Heisenberg antiferromagnet on a square lattice (cf. Fig. 

Clearly, the Hamiltonian then respects the full sym- 
metry of the lattice. Any of the couplings J c \, J C 2 and 
J aC 2 destroys this lattice isotropy, leading to a spatially 
anisotropic system for which cumulant ratios usually ex- 
hibit (small) deviations from their "isotropic" valuesi2£i2& 




FIG. 5: Binder cumulant U(80, 100), close to the boundary 
line of the antiferromagnetic phase, as a function of temper- 
ature. 



To determine the boundary of the spin-flop phase, we 
analysed the size dependence of the transverse sublattice 
magnetization, < Mj L>. We apply the criterion that the 
exponent g, see Eq. is 1/4 at the transition. Typical 
data are shown in Fig. demonstrating that the mag- 
netization decays much more rapidly with system size in 
the paramagnetic phase than in the spin-flop phase. To 
estimate the transition point, we determined the local 
slope (in a double logarithmic plot), 



dln<Mj > 



(9) 



from two consecutive system sizes, typically, L and L+20. 
Indeed, when crossing the phase boundary by fixing the 
temperature and lowering the field, g e g(L), for large L, 
tends to jump from 2, characterizing the decay in the 
disordered phase, to 1/4 at the transition to the spin- 
flop phase. Deeper in the spin-flop phase, g e s decreases 
slightly. 




FIG. 6: Logarithm of the transverse sublattice magnetization 
< M%y > versus the logarithm of the system size L at fixed 
temperature ksT = 0.6 and for various fields H z , close to the 
boundary between the spin-flop and the disordered phase. 



The Kosterlitz-Thouless character of the transition be- 
tween the spin-flop and the paramagnetic phases is also 
reflected in the thermal behavior of the specific heat C, 
which displays a non-critical maximum close to, but not 
exactly at the transition. Of course, from simulational 
data one cannot identify the expected essential singular- 
ity of C at the transition. 



IV. EFFECTS OF DEFECTS 

As discussed before^ experiments on La5CagCu2404i 
in a field H z along the easy axis provide no evidence for a 
sharp transition from the antiferromagnetic to the spin- 
flop phase. Instead, when fixing the temperature and in- 
creasing the field, the antiferromagnetic phase eventually 
becomes unstable against the disordered phase, and spin- 
flop structures seem to occur at higher fields only locally 
as indicated by a quite large, but non-critical maximum 
in the susceptibility)^ The reason for this experimentally 
observed behavior is not understood yet. Tentatively, one 
possible explanation invokes the holes or defects^ which 
may drive the transition and suppress the spin-flop phase. 

In the following, we shall explore this possibility by ex- 
tending the classical variant of the anisotropic Heisenberg 
model of Matsuda et al. by adding defects as described 
in Sec. [H] Actually, ten percent of the lattice sites will 
be occupied by these spinless, mobile defects, in accor- 
dance with the experiments^ We neglect the quantum 
nature of the holes and do not, e.g., include a kinetic en- 
ergy or "hopping" term in the Hamiltonian as one would 
normally expect in the case of a doped quantum antifer- 
romagnet. Of course, quantum effects may play an im- 
portant role for the phase behavior of the doped model. 
E.g., quantum fluctuations lead to a non-zero mobility of 
the holes even at T = where our spinless defects are 
static due to the absence of thermal fluctuations. Nev- 
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ertheless, our classical description is believed to provide 
some guidance to novel effects induced by the holes. 

Without external field (H z = 0) in the ground state 
(T = 0), the defects will form straight stripes perpen- 
dicular to the chains. Due to the next-nearest neighbor 
interactions, J C 2 and J aC 2, the stripes are bunched with 
two consecutive defects in a chain keeping the minimum 
distance of two lattice spacings with one spin in between 
the two defects. Such a bunching did not occur in the 
related Ising description, 13 where only nearest-neighbor 
couplings were assumed. In any event, the bunching may 
be suppressed, for example, by a pinning of the defects 
or repulsive interactions between the defects. The spins, 
in the ground state, are oriented along the z axis with an 
antiferromagnetic ordering from chain to chain, as in the 
case without defects. 

At low temperatures, the bunching dominates the typ- 
ical equilibrium configurations, as illustrated in Fig. E^. 
Upon increasing the temperature, the stripes tend to de- 
bunch, thereby gaining entropy, see Fig. [7|d. The de- 
bunching is reflected by a steep decrease in the density 
of defect pairs, i.e. consecutive defects in the same chain 
separated by merely one spin. The pronounced drop 
takes place in a rather narrow range of temperatures at 
roughly kgT « 0.55. However, the debunching seems 
to be a gradual, smooth process, without any thermal 
singularities. 

A phase transition occurs at UbT c « 0.7, i.e. at a signif- 
icantly lower temperature than in the absence of defects. 
At the transition, the defect stripes destabilize. As for 
the Ising model with mobile defects, the stripe instabil- 
ity may be inferred from the average minimal distance d a 
between each defect in chain to, at position (ld,m), and 
those in the next chain, at (l' d , m + 1), defined by 13-14 

d a = J2( min \ l d-l'd\)/N d , (10) 

dividing the sum by the number Nd of defects. This 
quantity increases rapidly at the transition. The transi- 
tion is also marked in the simulations by a pronounced 
peak in the specific heat and a drastic decrease in the 
sublattice magnetization, which is expected to vanish at 
T > T c in the thermodynamic limit. 

Applying an increasing external field H z > along 
the easy axis, the results of the ground state calcula- 
tion (for even numbers of at least four defects per chain) 
may be summarized as follows. First, one has to distin- 
guish two fields HP and hP . For < H z < hP one 
keeps the same antiferromagnetic structure with bunched 
defect stripes as in the case of vanishing field. Then, 
for hP < H z < hP, precisely one additional spin 
(pointing along the field direction) is inserted between 
two consecutive defects in every other chain, see Fig. 

(2) 

This configuration becomes unstable at H z — H z , and 
now an additional spin pointing in the z direction is in- 
serted between two consecutive defects in every chain 
(Fig. Eh). For simplicity and by analogy to the Ising 
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FIG. 7: Typical low temperature Monte Carlo equilibrium 
configurations showing (a) at fcsT = 0.3, the bunching of 
defect lines, and (b) at ksT = 0.6, the onset of debunching, 
both for systems of size L = 40. The open circles denote the 
mobile defects, while the arrows symbolize the z-components 
Sf m of the spins. 



caseji 3 - we refer to these two ground state configurations 
as "zig-zag" structures. The two fields HP and HP 
are readily found to be given by HP = —J ae 2 + J C 2/2 
and B~P = 2HP . Inserting the values of the interaction 
constants (see Sec. [HJ one obtains hP — 0.4305 and 

(2) 

H z ' — 0.861. Note that broader regions of "inserted" 
spins are, however, not favored energetically even when 
increasing the field. Instead, for larger fields, eventually 
a spin-flop transition occurs at H sf , followed by a ferro- 
magnetic structure at higher fields, as in the case without 
defects. 
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FIG. 8: Illustration of the ground state configurations in an 
external field H z > for a system with four defects per chain; 
for (a) H z 1] <H Z < H z 2 \ and (b) H z 2) < H z < H at . 



The zig-zag structures lead to a stepwise increase in 
the total magnetization, which gives rise to rather small 
maxima in the susceptibility x z at low temperatures, as 
depicted in Fig. [§] One observes a small, non-critical 
peak at H z « 0.9, which is the remnant of the transition 
at T = between the two zig-zag structures depicted 
in Fig. [S] Moreover, a very weak, non-critical maximum 
can be identified at H z w 0.5 (see inset). This peak is 
associated with the zig-zag structure of Fig.^t. At higher 
fields, a pronounced peak occurs signaling the transition 
to the spin-flop phase. 




FIG. 9: Susceptibility \ z as a function of the magnetic field 
H z at fixed temperature fesT = 0.3, for a system of size L = 
40. The inset shows the existence of a very weak maximum 
at H z » 0.5. 



We found, however, no evidence for a phase transition 
at finite temperatures associated with the small peaks in 
X z ■ Instead, upon increasing the field, straight stripes 
seem to transform gradually into zig-zag stripes, as for 



the Ising model with mobile defects^ The first peaks 
already vanishes at about kgT 0.3, and the position 
of the second small peak shifts to somewhat higher fields 
and gets less pronounced as the temperature is increased. 
It disappears at about ksT « 0.5, possibly due to the 
debunching. Obviously, the occurrence of zig-zag struc- 
tures cannot be identified with the phase transition in 
La5CagCu2404i observed well below the onset of spin- 
flop structures. 

Indeed, the anisotropic Heisenberg model with mobile 
defects still displays, at low temperatures, a sharp tran- 
sition from the antiferromagnetic phase, with straight or 
zig-zag defect stripes, to a spin-flop phase, as signaled by 
a delta-like peak in the susceptibility \ z ( see Fig- El- The 
topology of the phase diagram in the (T, H z ) plane seems 
to resemble that in the absence of defects, see Figs.Et and 
12)3. Actually, at the triple point (or critical end point) be- 
tween the antiferromagnetic, spin-flop and paramagnetic 
phases, located roughly at (fcsT = 0.5, H z = 1.7), the, 
presumably, non-critical debunching line seems to meet 
as well. However, we did not attempt to map the phase 
diagram accurately, because obviously the introduction 
of defects does not suffice to reconcile the experimental 
findings on La5CagCu2404i, showing no direct transi- 
tion from the antiferromagnetic to the spin-flop phase. 
In fact, the possible destruction of the spin-flop phase by 
the instability of the defect stripes tends to be hindered 
by the bunching of the stripes. Further investigations are 
desirable, but beyond the scope of the present study. 



V. SUMMARY 

We have analysed in detail a classical variant of a 
two-dimensional Heisenberg antiferromagnet with weak, 
uniaxial anisotropy proposed by Matsuda et al. to re- 
produce spin-wave dispersions measured in the magnet 
La5CagCu2404i. In particular, we determined the phase 
diagrams of the model applying fields parallel, H z , and 
perpendicular, H x , to the easy axis of the spin anisotropy. 

In the case of a transverse field H x (perpendicular to 
the easy axis), the transition from the antiferromagnetic 
phase to the paramagnetic phase belongs to the Ising uni- 
versality class. The phase diagram in the case of a field 
H z pointing along the easy axis consists of the antifer- 
romagnetic, the spin-flop, and the disordered (paramag- 
netic) phases. Extensive analyses have been performed 
to locate the phase boundaries, partly motivated by con- 
flicting analyses of related models. Indeed, our anal- 
ysis, studying especially the Binder cumulant and the 
transverse sublattice magnetization, allows one to locate 
reasonably well both the tricritical point on the phase 
boundary between the antiferromagnetic and the param- 
agnetic phases as well as the critical end point between 
these two phases and the spin-flop phase. Quantum ef- 
fects seem to play no essential role for the topology of 
the phase diagram, which is in qualitative disagreement 
with experimental observations on La5CagCu2404i. 
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We extended the classical variant of the model of Mat- 
suda et al. by including spinlcss mobile defects mimicking 
the holes in La5CagCu2404i, thereby following previous 
suggestions on a related Ising model. In the antiferro- 
magnetic phase, the defects, at low temperatures and low 
fields, arc found to form stripes as in the corresponding 
Ising case. However, due to next-nearest neighbor cou- 
plings, the stripes tend to bunch. The debunching, occur- 
ring at higher temperatures, seems to be non-critical, al- 
beit it takes place in a rather narrow temperature range. 
A phase transition at which the antiferromagnetic order 
is destroyed is driven by a destruction of the defect stripes 
loosing their coherency at the transition. 

The model with defects has also been studied in the 
presence of a field along the easy axis. There, a spin-flop 
phase is observed as well, separated from the antiferro- 
magnetic phase presumably by a transition of first order. 
Therefore, we conclude that adding the mobile defects 
is not sufficient to reconcile model properties with ex- 
perimental findings ruling out a direct transition from 
the antiferromagnetic to the spin-flop phase. Perhaps a 
destruction of the spin-flop phase may occur when the 
bunching is suppressed. 

However, when interpreting our findings for the model 
with defects one should keep in mind that our descrip- 



tion of the holes is a purely classical one. Quantum fluc- 
tuations reduce the clustering tendency of the holes and 
may also destroy the bunched structures that we find 
from our classical ground state analysis. Thus the role 
of quantum effects should certainly be investigated more 
carefully when comparing with actual experiments. 

In any event, the models display various interesting be- 
havior, and they may well contribute to arriving at a re- 
ally satisfying theoretical description of the intriguing ex- 
perimental observations on the La5CagCu2404i magnets. 
Moreover, the methods used in our study may be help- 
ful in analysing phase diagrams of other two-dimensional 
weakly anisotropic Heisenberg antiferromagnets. 
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